SUBROUTINE RDFVAR(CFILE,CVAR,POUTPUT)

USE YOMFORC1S, ONLY : JPSTPFC  , DTIMFC   ,RTSTFC   ,NSTPFC,DIMFORC
USE YOMLUN1S , ONLY : NULOUT   ,NULNAM
USE YOMDYN1S , ONLY : TSTEP
USE YOMCT01S , ONLY : NSTOP    ,NSTART
USE YOMGC1S  , ONLY : LMASK
USE YOMLOG1S , ONLY : NDIMCDF
USE YOMCST   , ONLY : RDAY 
USE YOMRIP   , ONLY : RTIMTR,RTIMST
USE YOMDPHY  , ONLY : NPOI,    NLON      ,NLAT     ,NLALO, NPOIP,NPOIPALL, NPOIALL, NPOIOFF
USE PARKIND1  ,ONLY : JPIM     ,JPRB, JPRD
USE NETCDF
USE NETCDF_UTILS, ONLY: NCERROR
USE YOMHOOK  ,ONLY : DR_HOOK, JPHOOK, LHOOK
USE MPL_MODULE
#ifdef DOC
! (C) Copyright 2000- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *RDFVAR  * - ROUTINE TO READ A SINGLE FIELD FROM A NETCDF FILE
!
! INTERFACE
! ---------
! CFILE		C	INPUT FILE
! CVAR		C	NETCDF VARIABLE NAME
! POUTPUT	R	ARRAY WITH OUTPUT VALUES
!
#endif
IMPLICIT NONE

#include "minmax.intfb.h"

CHARACTER(LEN=*),INTENT(IN) :: CFILE,CVAR
REAL(KIND=JPRB),DIMENSION(NPOI,JPSTPFC),INTENT(OUT) :: POUTPUT


!* LOCAL
CHARACTER(LEN=100) :: CNAME,CTUNITS
INTEGER(KIND=JPIM) ::  NCID,NIDTIM,NIDLON,NIDLAT,NITIM,NILAT,NILON,TVARID
INTEGER(KIND=JPIM) ::  VARID,VNDIMS,VDIMIDS(NF90_MAX_VAR_DIMS)
INTEGER(KIND=JPIM) ::  IYEAR,IMON,IDAY,IHOU,IMIN,ISEC,IFSS
INTEGER(KIND=JPIM) ::  IRDST,INSTFC,JT,NPTOT,NMX,NMY
INTEGER(KIND=JPIM) ::  ITUNITS(6)
REAL(KIND=JPRD)    :: ZFOR_START,ZFOR_REF,ZFOR_END,ZTSTRT,&
                      ZTIMEN,ZFCEN,ZZV,ZDELS

INTEGER ISTART4(4),ICOUNT4(4)
INTEGER ISTART3(3),ICOUNT3(3)
REAL(KIND=JPRB),ALLOCATABLE :: ZREAL3(:,:)
REAL(KIND=JPRB),ALLOCATABLE :: ZBUF(:)
INTEGER(KIND=JPIM) :: MYPROC, NPROC
INTEGER(KIND=JPIM) :: ISTP, IENP
REAL(KIND=JPHOOK)    :: ZHOOK_HANDLE

#include "fcttim.h"

IF (LHOOK) CALL DR_HOOK('RDFVAR',0,ZHOOK_HANDLE)

MYPROC = MPL_MYRANK()
NPROC  = MPL_NPROC()

ISTP  = NPOIOFF(MYPROC)+1
IENP = NPOIOFF(MYPROC+1)


POUTPUT(:,:)=0._JPRB
IF( MYPROC == 1 ) THEN
  !* Open forcing file 
  CALl NCERROR( NF90_OPEN(TRIM(CFILE),NF90_NOWRITE,NCID),'Opening '//TRIM(CFILE))
  WRITE(NULOUT,*)'NETCDF-FILE ',CFILE(1:len_trim(CFILE)),' OPENED ON UNIT ',NCID

  !! get variable ID
  CALL NCERROR( NF90_INQ_VARID(NCID, TRIM(CVAR), VARID),'getting varid '//TRIM(CVAR) )

  !* Find VARIABLE DIMENSIONS 
  CALL NCERROR( NF90_INQUIRE_VARIABLE(NCID,VARID,ndims=VNDIMS,dimids=VDIMIDS),'getting var info '//TRIM(CVAR) )
ENDIF

CALL MPL_BROADCAST(VARID,     KROOT=1,KTAG=100,CDSTRING='VARID')
CALL MPL_BROADCAST(VNDIMS,    KROOT=1,KTAG=200,CDSTRING='VNDIMS')
CALL MPL_BROADCAST(VDIMIDS(:),KROOT=1,KTAG=300,CDSTRING='VDIMIDS')


NIDTIM=VDIMIDS(VNDIMS)
NIDLAT=-1
IF( VNDIMS == 2 ) THEN
  NIDLON=VDIMIDS(1)
ELSEIF (VNDIMS == 3 .OR. VNDIMS == 4 ) THEN
  NIDLON=VDIMIDS(1)
  NIDLAT=VDIMIDS(2)
ELSE
  WRITE(NULOUT,*) 'WRONG NUMBER OF VARIABLE DIMENSIONS',TRIM(CFILE),TRIM(CVAR),VNDIMS
  CALL ABOR1('RDFVAR')
ENDIF
IF ( NIDLAT == -1 .AND. DIMFORC == 2 ) THEN
  WRITE(NULOUT,*) 'LAT NOT FOUND IN DIMENSION, BUT NAMELIST SAYS FORCING IS 2D'
  CALL ABOR1('RDFVAR')
ENDIF 

IF( MYPROC == 1 ) THEN
  !* get dimension lengths 
  CALL NCERROR( NF90_INQUIRE_DIMENSION(NCID,NIDTIM,len=NITIM),'getting time len '//TRIM(CVAR) )
  CALL NCERROR( NF90_INQUIRE_DIMENSION(NCID,NIDLON,len=NILON),'getting lon len '//TRIM(CVAR) )
  IF (NIDLAT /= -1 ) THEN
    CALL NCERROR( NF90_INQUIRE_DIMENSION(NCID,NIDLAT,len=NILAT),'getting lat len '//TRIM(CVAR) )
  ELSE
    NILAT = 0
  ENDIF

  IF(NILON.NE.NLON .OR. NILAT.NE.NLAT) THEN
    WRITE(*,*)'NLON OR NLAT NOT SPECIFIED CORRECTLY:'
    WRITE(*,*)'NLON IN NAMELIST AND FORCING: ',NLON,NILON
    WRITE(*,*)'NLAT IN NAMELIST AND FORCING: ',NLAT,NILAT
    CALL ABOR1('RDFVAR')
  ENDIF
ENDIF

CALL MPL_BROADCAST(NITIM,KROOT=1,KTAG=100,CDSTRING='NITIM')
CALL MPL_BROADCAST(NILON,KROOT=1,KTAG=200,CDSTRING='NILON')
CALL MPL_BROADCAST(NILAT,KROOT=1,KTAG=300,CDSTRING='NILAT')


IF( MYPROC == 1 ) THEN
!* Read time 
  CALL NCERROR( NF90_INQ_VARID(NCID, 'time', TVARID),'getting varid '//'TIME')
  CALL NCERROR( NF90_GET_ATT(NCID, TVARID, 'units', CTUNITS),'getting time units')
  PRINT*,CTUNITS
  IF ( CTUNITS(1:14) == "seconds since ") THEN
    READ(CTUNITS,'(14x,1i4,1x,1i2.2,1x,1i2.2,1x,1i2.2,1x,1i2.2,1x,1i2.2)')IYEAR,IMON,IDAY,IHOU,IMIN,ISEC
    ZDELS=1._JPRB
  ELSEIF ( CTUNITS(1:14) == "minutes since ") THEN
    READ(CTUNITS,'(14x,1i4,1x,1i2.2,1x,1i2.2,1x,1i2.2,1x,1i2.2,1x,1i2.2)')IYEAR,IMON,IDAY,IHOU,IMIN,ISEC
    ZDELS=60._JPRB
  ELSEIF ( CTUNITS(1:12) == "hours since ") THEN
    READ(CTUNITS,'(12x,1i4,1x,1i2.2,1x,1i2.2,1x,1i2.2,1x,1i2.2,1x,1i2.2)')IYEAR,IMON,IDAY,IHOU,IMIN,ISEC
    ZDELS=3600._JPRB
  ELSEIF ( CTUNITS(1:11) == "days since ") THEN
    READ(CTUNITS,'(11x,1i4,1x,1i2.2,1x,1i2.2,1x,1i2.2,1x,1i2.2,1x,1i2.2)')IYEAR,IMON,IDAY,IHOU,IMIN,ISEC
    ZDELS=86400._JPRB
  ELSE
    PRINT*, "Time units in netcdf not supported"
    PRINT*,CTUNITS
    CALL ABOR1('in rdfvar')
  ENDIF

  ITUNITS(1)=IYEAR
  ITUNITS(2)=IMON
  ITUNITS(3)=IDAY
  ITUNITS(4)=IHOU
  ITUNITS(5)=IMIN
  ITUNITS(6)=ISEC

  CALL NCERROR( NF90_GET_VAR(NCID,TVARID,ZFOR_START,start=(/1/)),'reading time start')
  CALL NCERROR( NF90_GET_VAR(NCID,TVARID,ZFOR_END,start=(/NITIM/)),'reading time end')

ENDIF

CALL MPL_BROADCAST(ZDELS,     KROOT=1,KTAG=100,CDSTRING='ZDELS')
CALL MPL_BROADCAST(ZFOR_START,KROOT=1,KTAG=200,CDSTRING='ZFOR_START')
CALL MPL_BROADCAST(ZFOR_END,  KROOT=1,KTAG=300,CDSTRING='ZFOR_END')
CALL MPL_BROADCAST(CTUNITS,   KROOT=1,KTAG=400,CDSTRING='CTUNITS')
CALL MPL_BROADCAST(ITUNITS(:),KROOT=1,KTAG=500,CDSTRING='ITUNITS')

IF( MYPROC > 1 ) THEN
  IYEAR= ITUNITS(1)
  IMON = ITUNITS(2)
  IDAY = ITUNITS(3)
  IHOU = ITUNITS(4)
  IMIN = ITUNITS(5)
  ISEC = ITUNITS(6)
ENDIF

!* forcing start/end 
IFSS=3600*IHOU+60*IMIN+ISEC
ZFOR_REF=RTIME(IYEAR,IMON,IDAY,IFSS)  !! forcing reference time 
ZFOR_START=ZFOR_REF+ZFOR_START*ZDELS  !! forcing start 
ZFOR_END=ZFOR_REF+ZFOR_END*ZDELS       !! forcing end 

!* check if forcing is available since model run start 
ZTSTRT=RTIMTR-0.5_JPRD*TSTEP  ! RTIMTR == ABSOLUTE TIME OF THE MODEL / Midpoint time of current time step (Julian seconds)
IF (ZFOR_START.GT.ZTSTRT) THEN
  WRITE(NULOUT,*) " STOP IN ROUTINE SUFCDF"
  WRITE(NULOUT,*) " RUN STARTS EARLIER THAN FORCING DATA"
  WRITE(NULOUT,*) " RUN STARTS AT ",ZTSTRT
  WRITE(NULOUT,*) " FORCING DATA STARTS ",ZFOR_START-ZTSTRT,&
     &     ' SECS LATER'
  WRITE(NULOUT,*) CTUNITS
  CALL ABOR1('in rdfvar')
ENDIF

!* find first forcing step needed to load:
! IRDST=NINT((ZTSTRT-ZFOR_START)/DTIMFC)+1
IRDST=FLOOR((ZTSTRT-ZFOR_START)/DTIMFC)+1 ! new
!RTSTFC=ZFOR_REF+DTIMFC*(IRDST-1)
RTSTFC=ZFOR_START+DTIMFC*(IRDST-1)
ZZV=MOD(NINT(RTSTFC-RTIME(IYEAR,IMON,IDAY,0)),NINT(RDAY))
IF ( ZZV  .NE. 0 ) THEN
   WRITE(*,*) '1st Forcing read not start at 00:00 UTC - this is problematic'
   WRITE(*,*) ' Please check this carefully to avoid interpolation problems!!' 
!    CALL ABOR1('in rdfvar')
ENDIF
!* find the number of steps to read 
! ZTIMEN=ZTSTRT+(NSTOP-NSTART+1)*TSTEP
ZTIMEN=ZTSTRT+(NSTOP-NSTART+1)*TSTEP 
! print*,'NSTOP,NSTART,IRDST,ZTIMEN',NSTOP,NSTART,IRDST,ZTIMEN
DO JT=IRDST,NITIM
  INSTFC=JT-IRDST+1
  ZFCEN=RTSTFC+DTIMFC*INSTFC
!   print*,'zz',JT,INSTFC,ZFCEN
  IF(ZFCEN.GT.ZTIMEN) EXIT
ENDDO
NSTPFC=INSTFC

!* check if forcing is available since model run start 
ZFCEN=RTSTFC+DTIMFC*(NSTPFC-1)
IF (ZFCEN+TSTEP.LT.ZTIMEN) THEN
  WRITE(NULOUT,*) " STOP IN ROUTINE SUFCDF"
  WRITE(NULOUT,*) " RUN ENDS LATER THAN FORCING DATA"
  WRITE(NULOUT,*) ' RUN BETWEEN ',ZTSTRT,' AND ',ZTIMEN
  WRITE(NULOUT,*) ' FORCING DATA BETWEEN ',RTSTFC,' AND ',ZFCEN,' Difference ',ABS(ZFCEN+TSTEP-ZTIMEN),' Time step model ',TSTEP
  CALL ABOR1('RDFVAR')
ENDIF

IF(NSTPFC.GT.JPSTPFC)THEN
  WRITE(*,*)'INCREASE JPSTPFC TO CONTAIN FORCING FIELDS'
  WRITE(*,*)'NSTPFC=',NSTPFC,'JPSTPFC=',JPSTPFC
  CALL ABOR1('RDFVAR')
ENDIF

WRITE(NULOUT,*) 'Forcing start-model start,Forcing end - model start,model end - model start (seconds since..):'
WRITE(NULOUT,*)RTSTFC-ZTSTRT,ZFCEN-ZTSTRT,ZTIMEN-ZTSTRT

WRITE(NULOUT,*)' FORCING READ BETWEEN FORCING TIME STEPS ',&
     &  IRDST,' AND ',IRDST+NSTPFC-1



!** ----- Read the data 
IF (DIMFORC == 2 ) THEN
  ISTART4(1:4) =(/1,1,1,IRDST/)
  ICOUNT4(1:4) =(/NILON,NILAT,1,NSTPFC/)
  ISTART3(1:3) =(/1,1,IRDST/)
  ICOUNT3(1:3) =(/NILON,NILAT,NSTPFC/)
  NPTOT=NILON*NILAT
  NMX=NLON
  NMY=NLAT
ELSE
  ISTART3(1:2) =(/1,IRDST/)
  ICOUNT3(1:2) =(/NILON,NSTPFC/)
  NPTOT=NILON
  NMX=NILON
  NMY=1
ENDIF

ALLOCATE (ZREAL3(NPTOT,NSTPFC))
ALLOCATE (ZBUF(NPOIALL))

IF( MYPROC == 1 ) THEN
  IF ( VNDIMS == 2 ) THEN
    CALL NCERROR( NF90_GET_VAR(NCID,VARID,ZREAL3,ISTART3(1:2),ICOUNT3(1:2)),'READING2 '//CVAR)
  ELSEIF ( VNDIMS == 3 ) THEN
    CALL NCERROR( NF90_GET_VAR(NCID,VARID,ZREAL3,ISTART3(1:3),ICOUNT3(1:3)),'READING3 '//CVAR)
  ELSEIF ( VNDIMS == 4 ) THEN
    CALL NCERROR( NF90_GET_VAR(NCID,VARID,ZREAL3,ISTART4(1:4),ICOUNT4(1:4)),'READING4 '//CVAR)
  ENDIF
  CALL NCERROR( NF90_CLOSE(NCID),'CLOSING FILE')
ENDIF

DO JT=1,NSTPFC
  CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3(:,JT),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDFVAR:POUTPUT")
  POUTPUT(1:NPOI,JT) = PACK(ZBUF(:),LMASK(ISTP:IENP))
ENDDO
IF( MYPROC == 1 ) THEN
  CALL MINMAX(CVAR,ZREAL3(:,1),NMX,NMY,LMASK,NULOUT)
ENDIF

WRITE(NULOUT,*) ' FORCING DATA ',TRIM(CVAR),' READ FOR ',NSTPFC,' FORCING STEPS'
IF (ALLOCATED(ZREAL3)) DEALLOCATE(ZREAL3)


IF (LHOOK) CALL DR_HOOK('RDFVAR',1,ZHOOK_HANDLE)

RETURN

END SUBROUTINE RDFVAR


